---
title: "Regressão Multinível"
author: "Patrícia Pain e Geovane Santos"
date: "`r Sys.Date()`"
output:
  html_document: default
  word_document: default
  pdf_document: default
editor_options:
  markdown:
    wrap: 72
---

# E quando o mundo não é linear?

# E quase sempre não é!

A modelagem linear cumpre bem o seu papel em determinados cenários, mas
em muitos casos ela sozinha não é suficiente para explicar um mundo tão
diverso. **E é ai que entra a modelagem multinível**! Um modelo que
prevê que o mundo é feito de grupos e que estes grupos exercem
influência sobre os dados. Este modelo considera que os dados podem ser
agrupados em até três níveis:

Onde:

-   Grupo 1 (t), podemos trabalhar com uma escala de tempo;
-   Grupo 2 (i), temos o indivíduo;
-   Grupo 3 (j) temos o agrupamento.

# Modelos Hierárquicos Lineares - Hierarchical Linear Modeling (HLM)

É uma forma complexa de regressão de mínimos quadrados ordinários (OLS)
usada para analisar a variância nas variáveis de resultado quando as
**variáveis preditoras estão em níveis hierárquicos variados**.

Destaque das limitações do OLS:

-   Ignora a estrutura hierárquica dos dados;

-   Supõe independência entre observações;

-   Não lida adequadamente com dados aninhados ou aninhados.

Razões para escolher a regressão multinível em vez do OLS: 

-   Modelagem de hierarquias naturais nos dados.

-   Controle da variabilidade intra e intergrupos.

-   Redução do viés e da inflação de erros padrão.

-   Melhor ajuste aos dados quando há dependência entre observações.

```         
Dados são separados em determinados grupos:
Por exemplo, dados podem ser separados por países, setor, ano 
```

As observações tendem a chegar mais próximo da realidade.

## Por que utilizar a modelagem multinível?

Modelo clássico de regressão aplicado a firmas de um mesmo país (sem que
hajam agrupamentos hierárquicos evidenciados pela literatura)

$y_{i} = \beta_{0} + \beta_{1}.X_{1i} + r_{i}$

Quando temos firmas de dois (ou mais) países diferentes, em que se
espera pela literatura uma relação entre as variáveis de interesse
distinta entre os países, duas (ou mais) equações distintas seriam
definidas

$y_{i1} = \beta_{01} + \beta_{11}.X_{1i1} + r_{i1}$

$y_{i2} = \beta_{02} + \beta_{12}.X_{1i2} + r_{i2}$

...

![](figura%20modelos%20clássicos.png)

A análise multinível é uma técnica mais adequada **para amostras com
estruturas agrupadas**, visto que reconhecem a presença de tais
hierarquias de dados, permitindo componentes residuais em cada nível na
hierarquia (Căpraru et al., 2018). De acordo com Fávero e Belfiore
(2017), isso **permite identificar e analisar as heterogeneidades
individuais e entre grupos a que pertencem estes indivíduos**, tornando
possível a **especificação de componentes aleatórios em cada nível** de
análise.

Os modelos multinível corrigem o fato de que as observações no mesmo
grupo não sejam independentes e, portanto, comparadas às modelagens
tradicionais, levam a estimativas não viesadas de erros padrão dos
parâmetros (Lima, Serra, Fávero, 2021). 

Ao comparar o modelo hierárquico multinível com outros modelos de
regressão ou com modelos de análise de covariância percebe-se que os
modelos hierárquicos multinível possuem como vantagem o fato de
considerarem a análise de dados hierarquicamente estruturados
(Goldstein, 2011; Hox, 2010).

Tal modelagem tem como benefício desagregar a variância da variável
dependente nos diversos níveis da análise, de maneira a **explicitar a
influência dos diversos níveis de análise** (Lima, Serra, Fávero,
2021). 

Conforme enfatizam Steenbergen e Jones (2002), Arceneaux e Nickerson
(2009) e Hair e Fávero (2019), se pesquisadores estiverem interessados
em testar se covariáveis em nível de grupo moderam os efeitos em nível
individual, os modelos multinível parecem ser a escolha mais apropriada.
E um teste de razão de verossimilhança pode ser elaborado para que se
verifique a adequação da estimação multinível à estrutura dos dados, bem
como propicia uma comparação com estimações oriundas, por exemplo, de
modelos tradicionais GLM. (Lima, Serra, Fávero, 2021). 

# Estrutura dos Dados

Nos modelos multinível, a estrutura dos dados nos permite combinar as
informações nos níveis individual e de grupo, isto é, entender o
**quanto e como a variação na relação entre variáveis se deve a
características individuais e/ou a efeitos 'contextuais'** (setores,
países, regiões, ...). As informações nos dois níveis se combinam
permitindo maior eficiência ao estimar coeficientes para grupos
específicos ou identificar efeitos no nível agregado (aninhamento).

Além de dados longitudinais, que agregam observações individuais no
tempo, dados com estrutura multinível envolvem dados aninhados. -\> Esse
tipo de dado tem um caráter hierárquico.

## Modelo Hierárquico Linear (HLM) de Dois Níveis

![](FIGURAS/Slide1.PNG)

Podemos pensar no desenho hierárquico equilibrado como em um painel
balanceado, já o desequilibrado se assemelha ao painel desbalanceado.

![](FIGURAS/Slide2.PNG)

### Modelo Hierárquico de Dois Níveis

$y_{ij} = \beta_{0j} + \beta_{1j}.X_{1ij} + r_{ij}$

Por termos dois níveis, $\beta_{0j}$ e $\beta_{1j}$ serão entendidos
como coeficientes aleatórios.

$\beta_{0j} = \gamma_{0} + u_{0j}$

$\beta_{1j} = \gamma_{1} + u_{1j}$

\*Gammas $\gamma's$ são os coeficientes de nível 2 e $u_{0j}$ e $u_{1j}$
são termos aleatórios com expectativa de 0.

**Modelo do Nível 1:**

$y_{ij} = \beta_{0j} + \beta_{1j}.X_{1ij} + \beta_{2j}.X_{2ij} + … + \beta_{Qj}.X_{Qij} + r_{ij}$

**Modelo do Nível 2:**

$\beta_{qj} = \gamma_{q0} + \gamma_{q1}.W_{1j} + \gamma_{q2}.W_{2j} + … + \gamma_{qS_{q}}.W_{S_{q}j} + u_{qj}$

## Modelo Hierárquico Linear (HLM) de Três Níveis

-   Aqui temos o mais utilizado nas pesquisas quando trabalhamos com
    dados em painel

![](FIGURAS/Slide3.PNG)

### As Medidas Repetidas

As variáveis de tempo, pois olhamos para empresas diferentes, de países
diferentes, em um mesmo ano. Diferente de quando tínhamos apenas dois
níveis, pois uma empresa não está em dois países ao mesmo tempo.

![](FIGURAS/Slide4.PNG)

### Modelo Hierárquico de Três Níveis

-   Grupo 1 (t) **Ano**, podemos trabalhar com uma escala de tempo;
-   Grupo 2 (i) **Firma**, temos o indivíduo;
-   Grupo 3 (j) **País**, temos o agrupamento.

Tínhamos:

<div>

$y_{ij} = \beta_{0j} + \beta_{1j}.X_{1ij} + r_{ij}$

Por termos dois níveis, $\beta_{0j}$ e $\beta_{1j}$ serão entendidos
como coeficientes aleatórios.

$\beta_{0j} = \gamma_{0} + u_{0j}$

$\beta_{1j} = \gamma_{1} + u_{1j}$

\*Gammas $\gamma's$ são os coeficientes de nível 2 e $u_{0j}$ e $u_{1j}$
são termos aleatórios com expectativa de 0.

**Modelo do Nível 1:**

$y_{ij} = \beta_{0j} + \beta_{1j}.X_{1ij} + \beta_{2j}.X_{2ij} + … + \beta_{Qj}.X_{Qij} + r_{ij}$

**Modelo do Nível 2:**

$\beta_{qj} = \gamma_{q0} + \gamma_{q1}.W_{1j} + \gamma_{q2}.W_{2j} + … + \gamma_{qS_{q}}.W_{S_{q}j} + u_{qj}$

</div>

E agora:

**Modelo do Nível 1:**

$y_{tij} = \pi_{0ij} + \pi_{1ij}.ANO_{tij} + e_{yij}$

$\pi_{0ij}$: valor esperado da variável dependente da firma ij no ano 1;

$\pi_{1ij}$: taxa de crescimento da variável dependente da firma ij;

$e_{tij}$: variância da firma ao longo do tempo.

**Modelo do Nível 2:**

$\pi_{pij} = \beta_{p0j} + \beta_{p1j}.X_{1ij} + \beta_{p2j}.X_{2ij} + … + \beta_{pQ_{p}j}.X_{Q_{p}ij} + r_{pij}$

<div>

$\pi_{1ij}$: taxa de crescimento da variável dependente da firma ij;

</div>

$\beta$ (q = 0, 1, ..., $Q_{p}$): coeficientes do nível 2;

$X$: vetor das variáveis preditoras do nível 2;

$r_{pij}$: efeito aletório do nível 2.

**Modelo do Nível 3:**

$\beta_{pqj} = \gamma_{pq0} + \gamma_{pq1}.W_{1j} + \gamma_{pq2}.W_{2j} + … + \gamma_{pqS_{pq}}.W_{S_{pq}j} + u_{pqj}$

<div>

$\beta$ (q = 0, 1, ..., $Q_{p}$): coeficientes do nível 2;

</div>

$\gamma$ (s = 0, 1, ..., $S_{pq}$): coeficientes do nível 3;

$W$: vetor das variáveis preditoras do nível 3;

$u_{pqj}$: efeito aletório do nível 3.

# A Importância da Teoria

| Nível 1 | Nível 2 | Nível 3 |
|---------|---------|---------|
| Ticker  | Setor   |         |
| Ano     | Ticker  | Setor   |
| Ticker  | País    |         |
| Ticker  | Setor   | País    |
| Ano     | Ticker  | País    |
| ?       | ?       | ?       |

Os dados também podem ser não aninhados, como ocorre em modelos que
inferem a posição ideológica de deputados por meio de seu padrão de
votação nominal e onde temos votos agregados no nível dos deputados e do
projeto de lei.

Podemos pensar nos modelos multinível como uma regressão onde incluímos
indicadores para cada grupo.

## Parâmetros

# Como Operacionalizar no R?

## Instalar e carregar pacotes

-   plotly

-   tidyverse

-   reshape2

-   knitr

-   kableExtra

-   rgl

-   car

-   nlme

-   lmtest

-   fastDummies

-   msm

-   lmeInfo

-   jtools

-   data.table

-   olsrr

```{r pacotes, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
pacotes <- c("basictabler","car","caret","conflicted","correlation","cowplot", "data.table","dplyr","faraway","fastDummies","flextable", "foreign","ggrepel","ggtree","ggplot2","ggpubr","graphics", "grid", "gridExtra","gtsummary","Hmisc","jsmodule", "jtools","knitr","kableExtra","lmtest", "lmeInfo", "lubridate","magick", "MASS","MatchIt","mfx","minqa","mgcv", "msm", "nlme","nnet", "olsrr", "palmerpenguins","performance","pglm",  "pgls","plm","plotly","pROC","pscl","psych","readr","readxl","ReporteRs","reshape2", "rgl", "ROCR","scales","sjlabelled","stargazer","stats","tidyr","tidyverse","tinytex",  "tseries","vdr","visreg","xlsx","xtable","writexl")
```

```{r carregando, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
options(repos = "https://cran.rstudio.com/")
if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
  instalador <- pacotes[!pacotes %in% installed.packages()]
  for(i in 1:length(instalador)) {
    install.packages(instalador, dependencies = T)
    break()}
  sapply(pacotes, require, character = T)
} else {
  sapply(pacotes, require, character = T)
}
```

## Função para significância dos níveis

A função `stderr_nlme`, de autoria dos professores Fávero e Rafael Sousa
ajuda a identificar a significância estatística dos dados de cada nível
e de certa forma a calcular o quanto do comportamento dos dados é
explicado por cada nível.

```{r}
stderr_nlme <- function(model){
  if(base::class(model) != "lme"){
    base::message("Use a lme object model from nlme package")
    stop()}
  resume <- base::summary(model)
  if(base::length(base::names(model$groups))==1){
    m.type <- "HLM2"
  } else if(base::length(base::names(model$groups))==2){
    m.type <- "HLM3"
  }
  if(m.type == "HLM2"){
    vcov_matrix <- model$apVar
    logs_sd_re <- base::attr(vcov_matrix,"Pars")
    if(base::length(logs_sd_re)==2){
      stderr_tau00 <- msm::deltamethod(~exp(x1)^2,logs_sd_re,vcov_matrix)
      stderr_sigma <- msm::deltamethod(~exp(x2)^2,logs_sd_re,vcov_matrix)
      results <- base::data.frame(`RE Components`=base::c("Var(v0j)","Var(e)"),
                                  `Variance Estimatives`= base::c(base::exp(logs_sd_re)[[1]]^2,
                                                                  base::exp(logs_sd_re[[2]])^2),
                                  `Std Err.`=base::c(stderr_tau00,
                                                     stderr_sigma),
                                  z=base::c(base::exp(logs_sd_re)[[1]]^2/stderr_tau00,
                                            base::exp(logs_sd_re[[2]])^2/stderr_sigma),
                                  `p-value`=base::round(stats::pnorm(q=base::c(base::exp(logs_sd_re)[[1]]^2/stderr_tau00,
                                                                               base::exp(logs_sd_re[[2]])^2/stderr_sigma),
                                                                     lower.tail=F)*2,3))
      return(results)
    }
    else{
      stderr_tau00 <- msm::deltamethod(~exp(x1)^2,logs_sd_re,vcov_matrix)
      stderr_tau01 <- msm::deltamethod(~exp(x2)^2,logs_sd_re,vcov_matrix)
      stderr_sigma <- msm::deltamethod(~exp(x4)^2,logs_sd_re,vcov_matrix)
      results <- base::data.frame(Components=base::c("Var(v0j)","Var(v1j)","Var(e)"),
                                  `Variance Estimatives`= base::c(base::exp(logs_sd_re)[[1]]^2,
                                                       base::exp(logs_sd_re[[2]])^2,
                                                       base::exp(logs_sd_re[[4]])^2),
                                  `Std Err.`=base::c(stderr_tau00,
                                                  stderr_tau01,
                                                  stderr_sigma),
                                  z=base::c(base::exp(logs_sd_re)[[1]]^2/stderr_tau00,
                                            base::exp(logs_sd_re[[2]])^2/stderr_tau01,
                                            base::exp(logs_sd_re[[4]])^2/stderr_sigma),
                                  `p-value`=base::round(stats::pnorm(q=base::c(base::exp(logs_sd_re)[[1]]^2/stderr_tau00,
                                                                               base::exp(logs_sd_re[[2]])^2/stderr_tau01,
                                                                               base::exp(logs_sd_re[[4]])^2/stderr_sigma),
                                                                     lower.tail=F)*2,3))
      return(results)
    }
  }
  if(m.type == "HLM3"){
    vcov_matrix <- model$apVar
    logs_sd_re <-  base::attr(vcov_matrix,"Pars")
    if(base::length(logs_sd_re) == 3){
      stderr_tau_r000 <- msm::deltamethod(~exp(x1)^2,logs_sd_re,vcov_matrix)
      stderr_tau_u000 <- msm::deltamethod(~exp(x2)^2,logs_sd_re,vcov_matrix)
      stderr_sigma <- msm::deltamethod(~exp(x3)^2,logs_sd_re,vcov_matrix)
      results <- base::data.frame(Components=base::c("Var(t00k)","Var(v0jk)","Var(e)"),
                                  Estimatives=base::c(base::exp(logs_sd_re)[[2]]^2,
                                                      base::exp(logs_sd_re)[[1]]^2,
                                                      base::exp(logs_sd_re)[[3]]^2),
                                  Std_Err=base::c(stderr_tau_u000,
                                                  stderr_tau_r000,
                                                  stderr_sigma),
                                  z=base::c(base::exp(logs_sd_re)[[2]]^2/stderr_tau_u000,
                                            base::exp(logs_sd_re)[[1]]^2/stderr_tau_r000,
                                            base::exp(logs_sd_re)[[3]]^2/stderr_sigma),
                                  `p-value`=base::round(stats::pnorm(q=base::c(base::exp(logs_sd_re)[[2]]^2/stderr_tau_u000,
                                                                               base::exp(logs_sd_re)[[1]]^2/stderr_tau_r000,
                                                                               base::exp(logs_sd_re)[[3]]^2/stderr_sigma),
                                                                     lower.tail=F)*2,3))
      return(results)
    } 
    else{
      stderr_tau_r000 <- msm::deltamethod(~exp(x1)^2,logs_sd_re,vcov_matrix)
      stderr_tau_r100 <- msm::deltamethod(~exp(x2)^2,logs_sd_re,vcov_matrix)
      stderr_tau_u000 <- msm::deltamethod(~exp(x4)^2,logs_sd_re,vcov_matrix)
      stderr_tau_u100 <- msm::deltamethod(~exp(x5)^2,logs_sd_re,vcov_matrix)
      stderr_sigma <- msm::deltamethod(~exp(x7)^2,logs_sd_re,vcov_matrix)
      results <- base::data.frame(`RE_Components`=base::c("Var(t00k)","Var(t10k)",
                                                          "Var(v0jk)","Var(v1jk)",
                                                          "Var(e)"),
                                  `Variance Estimatives`=base::c(base::exp(logs_sd_re)[[4]]^2,
                                                                 base::exp(logs_sd_re)[[5]]^2,
                                                                 base::exp(logs_sd_re)[[1]]^2,
                                                                 base::exp(logs_sd_re)[[2]]^2,
                                                                 base::exp(logs_sd_re)[[7]]^2),
                                  `Std Err.`=base::c(stderr_tau_u000,
                                                     stderr_tau_u100,
                                                     stderr_tau_r000,
                                                     stderr_tau_r100,
                                                     stderr_sigma),
                                  z=base::c(base::exp(logs_sd_re)[[4]]^2/stderr_tau_u000,
                                            base::exp(logs_sd_re)[[5]]^2/stderr_tau_u100,
                                            base::exp(logs_sd_re)[[1]]^2/stderr_tau_r000,
                                            base::exp(logs_sd_re)[[2]]^2/stderr_tau_r100,
                                            base::exp(logs_sd_re)[[7]]^2/stderr_sigma),
                                  `p-value`=base::round(stats::pnorm(q=base::c(base::exp(logs_sd_re)[[4]]^2/stderr_tau_u000,
                                                                               base::exp(logs_sd_re)[[5]]^2/stderr_tau_u100,
                                                                               base::exp(logs_sd_re)[[1]]^2/stderr_tau_r000,
                                                                               base::exp(logs_sd_re)[[2]]^2/stderr_tau_r100,
                                                                               base::exp(logs_sd_re)[[7]]^2/stderr_sigma),
                                                                     lower.tail=F)*2,3))
      return(results)
    }
  }
}
```

## Carregando a Base

Base para análise e Base teste para comparação entre previsto e real

```{r}
load("C:/Users/Usuário/OneDrive/Oficina multinivel/Environment inicial.RData")
```

### Análise exploratória

Ao visualizarmos a função de densidade de probabilidade da função
dependente (Qualidade da Auditoria) observa-se distribuições diferentes
dependendo do país, mostrando que temos uma dispersão grande dependendo
do país.

```{r}
Base_HLM %>% 
  group_by(Country) %>% 
  mutate(linhas = 1:n()) %>% 
  mutate(x = unlist(density(ACCY_w1, n = max(linhas))["x"]),
         y = unlist(density(ACCY_w1, n = max(linhas))["y"])) %>%
  mutate_if(is.numeric, scale) %>%
  ggplot() +
  geom_area(aes(x = x, y = y, group = Country, fill = Country), color = "black", alpha = 0.4) +
  geom_histogram(aes(x = ACCY_w1, 
                     y = ..density.., 
                     fill = Country), 
                 color = "black", 
                 position = 'identity', 
                 alpha = 0.3) +
  facet_wrap(~ Country) +
  scale_fill_viridis_d() +
  scale_color_viridis_d() +
  theme_classic()
```

## Será o HLM o modelo mais adequado?

Para a modelagem multinível utiliza-se a técnica **step-up** onde o
modelo é otimizado a cada execução até obter-se sua versão final.

Para verificar o modelo vamos utilizar o **BIC (bayesian information
criterion)** critério que avalia o quão perto o modelo está dos dados
reais considerando também as variáveis explicativas.

### Modelo nulo

O primeiro modelo a ser criado é o **modelo nulo**, que não considera
nenhuma variável preditora, apenas a ACCY e o agrupamento por país.

No R a execução desta fórmula é feita pelo comando `lme` do pacote
`nlme`

*Olhando para o nível País*

```{r}
hlm_nulo = 
  lme(fixed = ACCY_w1 ~1, 
      random = ~1 | Country, 
      data = Base_HLM,
      method = "REML")
```

O BIC do modelo base será utilizado mais adiante para compararmos o
ganho dos próximos modelos. Quando menor ele for, melhor!

```{r}
BIC(hlm_nulo)
```

| Modelo | BIC      |
|--------|----------|
| Nulo   | 77044.39 |

```{r}
stderr_nlme(hlm_nulo)
```

O p-value de nossa variável ACCY (v0j) é estatisticamente significativa
(menor que 0,05), indicando que o modelo multinível é aplicável à este
caso. Se ela fosse superior a 0,05 o modelo OLS seria o mais indicado.

Utilizaremos agora a base de testes para executar o modelo e comparar se
os resultados previstos estão próximos do resultado real.

\*A base de testes (teste) é uma base idêntica.

```{r}
teste$hlm_nulo_fitted = predict(hlm_nulo, teste)

teste %>%
  mutate_if(is.numeric, scale) %>%
  ggplot() +
  geom_smooth(aes(x = ACCY_w1, y = hlm_nulo_fitted, color = "1 - Modelo Nulo"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +
  geom_smooth(aes(x = ACCY_w1, y = ACCY_w1), method = "lm", 
              color = "gray44", size = 1.05,
              linetype = "longdash") +
  labs(x = "ACCY real", y = "ACCY Predito") +
  theme_bw()
```

Observe que a linha rosa (modelo nulo) chega a cruzar a linha pontilhada
(resultado perfeito), mas apresenta considerável diferença.

### Modelo com interceptos aleatórios

O próximo passo, é o **modelo com interceptos aleatórios** que considera
apenas efeito que o país (intercepto) exerce sobre a ACCY:

```{r}
hlm_intercept_aleatorios = 
  lme(fixed = ACCY_w1 ~  SIZE_w1  + Emerg,
      random = ~1 | Country,
      data = Base_HLM,
      method = "REML")
```

Observação do BIC do modelo. Se o indicador BIC for reduzido em relação
ao modelo nulo indica que a variável Tamanho possui grande efeito sobre
a qualidade da auditoria (ACCY).

```{r}
BIC(hlm_intercept_aleatorios)
```

| Modelo                 | BIC      |
|------------------------|----------|
| Nulo                   | 77044.39 |
| Interceptos aleatórios | 76701.39 |

```{r}
stderr_nlme(hlm_intercept_aleatorios)
```

O *p-value* do nível v0j abaixo de 0,05 indica que o modelo é
estatisticamente significativo.

Vamos verificar agora como os modelos se aproximam do valor ideal:

```{r}
teste$hlm_intercept_aleatorios_fitted = predict(hlm_intercept_aleatorios, teste)


teste %>%
  mutate_if(is.numeric, scale) %>%
  ggplot() +
  geom_smooth(aes(x = ACCY_w1, y = hlm_nulo_fitted, color = "1 - Modelo Nulo"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +
  geom_smooth(aes(x = ACCY_w1, y= hlm_intercept_aleatorios_fitted, color = "2 - Interceptos Aletaórios"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +
  geom_smooth(aes(x = ACCY_w1, y = ACCY_w1), method = "lm", 
              color = "gray44", size = 1.05,
              linetype = "longdash") +
  labs(x = "ACCY real", y = "ACCY Predito") +
  theme_bw()
```

Observe que a linha do modelo com interceptos aleatórios se aproxima
mais da reta ideal do que a linha do modelo nulo.

### Modelo com interceptos e inclinações aleatórias

O próximo passo é verificarmos o efeito que o país exerce sobre o
Tamanho e a ACCY.

```{r}
hlm_intercept_inclin = 
  lme(fixed = ACCY_w1 ~ SIZE_w1  + Emerg, 
      random = ~ SIZE_w1  | Country, 
      data = Base_HLM,
      method = "REML")
```

```{r}
BIC(hlm_intercept_inclin)
```

Nesta terceira execução observa-se uma nova queda no indicador BIC,
mostrando que o volume de erros diminuiu em comparação com os modelos
anteriores

| Modelo                               | BIC      |
|--------------------------------------|----------|
| Nulo                                 | 77044.39 |
| Interceptos aleatórios               | 76701.39 |
| Interceptos e Inclinações aleatórias | 75948.27 |

```{r}
stderr_nlme(hlm_intercept_inclin)
```

O *p-value* do nível v0j (nível firma) e v1j se seguirem abaixo de 0,05
indicam que o modelo continua estatisticamente significativo.

Vamos comparar a performance dos modelos

```{r}
teste$hlm_intercept_inclin_fitted = predict(hlm_intercept_inclin, teste)


teste %>%
  mutate_if(is.numeric, scale) %>%
  ggplot() +
  geom_smooth(aes(x = ACCY_w1, y = hlm_nulo_fitted, color = "1 - Modelo Nulo"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +
  geom_smooth(aes(x = ACCY_w1, y= hlm_intercept_aleatorios_fitted, color = "2 - Interceptos Aletaórios"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +
  geom_smooth(aes(x = ACCY_w1, y= hlm_intercept_inclin_fitted, color = "3 - Interceptos e Inclinações aleatórias"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +  
  geom_smooth(aes(x = ACCY_w1, y = ACCY_w1), method = "lm", 
              color = "gray44", size = 1.05,
              linetype = "longdash") +
  labs(x = "ACCY real", y = "ACCY Predito") +
  theme_bw()
```

A linha do modelo com interceptos e inclinações aleatórias se aproxima
mais ainda da reta ideal do que as linhas dos demais modelos.

### Modelo Final

Agora que verificamos que o Tamanho tem efeito sobre a ACCY e que o
agrupamento por países também exerce efeito sobre as demais informações,
vamos verificar se a variável no nível de país também impacta no
resultado. Teremos finalmente a fórmula completa:

```{r}
hlm_final = 
  lme(fixed = ACCY_w1 ~ SIZE_w1  + Law + SIZE_w1:Law,
      random = ~ SIZE_w1 | Country,
      data = Base_HLM,
      method = "REML")
```

```{r}
BIC(hlm_final)
```

Se tiver uma nova redução no BIC, indica que o houve melhora no modelo
ao incluir o efeito do Tamanho no sistema legal e na ACCY.

| Modelo                               | BIC      |
|--------------------------------------|----------|
| Nulo                                 | 77044.39 |
| Interceptos aleatórios               | 76701.39 |
| Interceptos e Inclinações aleatórias | 75948.27 |
| Final                                | 75958.33 |

```{r}
stderr_nlme(hlm_final)
```

O *p-value* do nível v0j (nível da firma) e v1j seguem abaixo de 0,05 ou
seja, modelo significativo estatisticamente.

Vamos verificar visualmente como essa nova reta se porta com a base de
testes:

```{r}
teste$hlm_final_fitted = predict(hlm_final, teste)


teste %>%
  mutate_if(is.numeric, scale) %>%
  ggplot() +
  geom_smooth(aes(x = ACCY_w1, y = hlm_nulo_fitted, color = "1 - Modelo Nulo"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +
  geom_smooth(aes(x = ACCY_w1, y= hlm_intercept_aleatorios_fitted, color = "2 - Interceptos Aletaórios"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +
  geom_smooth(aes(x = ACCY_w1, y= hlm_intercept_inclin_fitted, color = "3 - Interceptos e Inclinações aleatórias"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +  
  geom_smooth(aes(x = ACCY_w1, y= hlm_final_fitted, color = "4 - Modelo Final"),
              method = "lm", se = F, formula = y ~ splines::bs(x, df = 5),
              size = 1.5) +  

  geom_smooth(aes(x = ACCY_w1, y = ACCY_w1), method = "lm", 
              color = "gray44", size = 1.05,
              linetype = "longdash") +
  labs(x = "ACCY real", y = "ACCY Predito") +
  theme_bw()
```

Ao compararmos as retas de regressão observamos que o modelo final
sobrepõe quase que perfeitamente a reta de interceptos e inclinações
aleatórias, ao considerar o efeito de todas as variáveis disponíveis.

## HLM x OLS

### Dois Níveis

Rodando o modelo nulo

```{r}
install.packages("glmmTMB")
library(glmmTMB)
```

```{r}
hlm2_nullmodel <- glmmTMB(formula = ACCY_w1 ~ 1 + (1 | Ticker), family = gaussian, data = Base_HLM, REML = TRUE)

summary(hlm2_nullmodel)
```

Olhando os interceptos para firmas

```{r}
ranef(object = hlm2_nullmodel)
```

Rodando o OLS nulo para comparar com o HLM

```{r}
ols_nullmodel <- lm(formula = ACCY_w1 ~ 1, data = Base_HLM)

summary(ols_nullmodel)
```

Comparando os modelos, se `hlm2_nullmodel` for significativo a 1%, então
HLM é melhor que OLS nesse caso

```{r}
lrtest(ols_nullmodel, hlm2_nullmodel)
```

#### Rodando o modelo final

Dentro dos parênteses, declaramos **variável de interesse \| variável de
nível**, determinando que o R calcule os interceptos aleatórios causados
pela variável de interesse devido aos diferentes participantes dos
níveis (firmas).

```{r}
hlm2_finalmodel <- glmmTMB(formula = ACCY_w1 ~ SIZE_w1 + Emerg + Law + CFO_w1 + AU + (SIZE_w1 | Ticker), 
                           family = gaussian, 
                           data = Base_HLM, 
                           REML = TRUE)
summary(hlm2_finalmodel)
```

### Três Níveis

\>Sempre o maior nível primeiro

Rodando o modelo nulo

```{r}
hlm3_nullmodel <- glmmTMB(formula = ACCY_w1 ~ 1 + (1 | Country) + (1 | Ticker), family = gaussian, data = Base_HLM, REML = TRUE)

summary(hlm3_nullmodel)
```

Olhando os interceptos para firmas e países

```{r}
ranef(object = hlm3_nullmodel)
```

Rodando o OLS nulo para comparar com o HLM

```{r}
#O mesmo para 2 ou 3 níveis
ols_nullmodel <- lm(formula = ACCY_w1 ~ 1, data = Base_HLM)
```

Comparando os modelos, se `hlm3_nullmodel` for significativo a 1%, então
HLM é melhor que OLS nesse caso

```{r}
lrtest(ols_nullmodel, hlm3_nullmodel)
```

#### Rodando o modelo final

Dentro dos parênteses, declaramos **variável de interesse \| variável de
nível**, determinando que R calcule os interceptos aleatórios causados
pela variável de interesse devido aos diferentes participantes dos
níveis (firmas e países).

```{r}
hlm3_finalmodel <- glmmTMB(formula = ACCY_w1 ~ SIZE_w1 + Emerg + Law + CFO_w1 + AU + (SIZE_w1 | Country) + (SIZE_w1 | Ticker), 
                           family = gaussian, 
                           data = Base_HLM, 
                           REML = TRUE)
summary(hlm3_finalmodel)
```
